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SECTION  I 

MICROTHERMAL  SPECTRAL  EXAMINATION 


This  report  is  concerned  with  the  estimation  of  power  spectra  of 
microthermal  fluctuations.  It  represents  work  directed  at  verifying 
spectra  being  routinely  obtained  from  microthermal  data  recorded 
using  high-flying  aircraft.  The  material  within  this  report  is  not 
intended  to  represent  new  and  unique  research  in  the  field  of  digital 
signal  processing,  but  rather  to  document  and  discuss  the  various 
well  accepted  techniques. 

The  microthermal  fluctuations  are  important  because  they  give  in- 
formation on  spatial  refractive  index  fluctuations  which  is  necessary 
to  the  calculation  of  the  behavior  of  light  beams  propagating  to  and 
from  high-flying  aircraft.  During  1974  and  1975  many  such  spectra  were 
obtained  at  Kirtland  Air  Force  Base  to  provide  a heretofore  unavailable 
data  base  for  such  calculations.  The  calculations  presented  here  were 
performed  as  a routine  check  on  these  spectra  and  the  procedures  used 
to  obtain  them. 

The  microthermal  fluctuation  data  are  obtained  using  a platinum- 
wire  microthermal  sensor  mounted  near  the  front  or  top  on  an  aircraft 
fuselage  and  sufficiently  removed  from  the  surface  to  be  outside  the 
aircraft  boundary  layer.  The  data  are  recorded  on  FM  analog  tape  in 
sections  roughly  two  minutes  in  duration.  Typically  six  to  twelve 
such  runs  might  be  recorded  in  one  flight.  The  data  are  then  possibly 
filtered,  digitized,  and  computer  processed  to  obtain  the  temporal 
spectra.  The  spectra  obtained  here  check  on  the  data  processing  portion 
of  the  effort,  repeating  the  digitization  and  computer  processing. 

Section  II  of  this  report  contains  a discussion  of  the  preliminary 
data  management  including  digitization  rates,  scaling  and  the  removal 
of  some  spikes  that  were  found  in  the  analog  record.  Sections  III  and 
IV  contain  the  details  of  the  procedures  for  obtaining  the  spectra. 

The  spectra  were  divided  into  frequency  regions.  Section  III  contains 
the  details  of  the  basic  processing  procedure.  These  were  used  di- 
rectly in  the  highest  frequency  region.  The  basic  procedure  includes 
multiplying  by  data  windows,  pre-whitening/post-darkening,  ensemble 
averaging,  and  the  calculation  of  confidence  intervals.  In  Section  IV 
the  additional  procedures  of  smoothing  and  decimation,  necessary  for 
extending  the  spectrum  to  lower  frequency  ranges,  are  presented.  The 
results  of  the  combined  high  and  low  frequency  procedures  are  presented 
in  Section  V.  Section  V also  contains  a comparison  with  the  spectrum 
obtained  previously  at  Kirtland  Air  Force  Base.  Section  VI  gives  some 
thoughts  on  procedures  suggested  to  improve  spectral  resolution  and 
rel iabil ity. 
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In  this  section  we  present  a discussion  of  the  sampling  rate  and 
a discussion  of  one  other  problem,  the  occurence  of  spikes  in  the 
primary  data  tape. 

To  be  thorough  we  elected  to  determine  our  own  digitization  rate. 
The  full  bandwidth  of  the  electronics  was  20  KHz.  To  achieve  at  least 
this  sampling  frequency  with  the  available  500  Hz  A-D  converter,  we 
played-back  the  tape  at  a speed  reduced  from  30  IPS  by  a factor  of  64, 
i.e.,  at  a tape  speed  of  15/32  IPS.  At  this  tape  speed  the  Nyquist 
rate  corresponding  to  our  system  bandwidth  is  312.5  Hz. 

Calibration  of  the  data  was  accomplished  by  means  of  a -.15°C 
100  Hz  square  wave  recorded  on  the  leading  edge  of  the  tape.  Several 
cycles  of  this  reference  square  wave  were  digitized  along  with  approxi- 
mately 2 minutes  of  data.  In  addition  to  calibrating  the  magnitude  of 
the  temperature  fluctuations,  inspection  of  the  digitized  version  of 
the  square  wave  enabled  the  accurate  calculation  of  the  sampling  rate. 
This  precise  value  of  the  equivalent  sampling  rate  was  30.65  KHz.  In- 
spection of  the  digitized  microthermal  data  revealed  several  periods 
with  severe  spiking.  Such  spiking  is  exemplified  by  Figure  1,  which 
represents  33.5  msec,  of  data.  (The  straight  line  represents  the 
average  temperature.)  The  spiking  in  general  consisted  of  from  two 
to  five  very  narrow  (^130  ysec)  pulses  spaced  approximately  8 msec, 
apart.  The  inter-group  spacing  appeared  to  be  random.  Since  in- 
clusion of  these  spikes  in  our  spectral  calculations  would  decrease 
the  reliability  of  our  results,  they  were  eliminated  by  a simple  linear 
interpolation  between  data  points  on  either  side  of  the  spikes.  The 
same  time  span  of  data  shown  in  Figure  1 is  again  shown  in  Figure  2 
with  the  spikes  removed  by  the  interpolation  scheme. 

This  concludes  our  discussion  of  the  preliminary  processing  of  the 
microthermal  data  including  digitization,  scaling  and  spike  removal. 

The  next  section  outlines  the  techniques  employed  in  the  hi gh-frequency 
analysis. 
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Figure  1.  Time  trace  of  raw  microthermal  data 


Figure  2.  Time  trace  of  microthermal  data 
with  spikes  removed. 


SECTION  III 
HIGH  FREQUENCY  RANGE 

We  first  proceed  with  a discussion  of  some  classical  techniques 
employed  in  the  digital  estimation  of  power  spectra. 

A.  Power  Spectrum  Using  the  Fast  Fourier  Transform 

Here  we  state  the  equations  used  in  the  Fast  Fourier  Transform 
(FFT)  and  in  the  calculation  of  the  discrete  power  spectrum  and  show 
that  they  approximate  the  continuous  spectral  integrals  under  consid- 
eration. We  assume  that  we  have  many  groups  of  N sampled  data  values, 
sampled  at  intervals  of  tQ. 

The  equations  implemented  by  the  FFT  algorithm  are 
N i I1  (K-1)(L-1) 

X (L ) = l C(K)  e N (la) 


1 N -i  jr  (K-D(L-l) 

C(K)  = J jj  X(L)e 


This  pair  of  equations  can  be  shown  to  approximate  the  continuous 
Fourier  transform  pair 


57  g(u))e1wtdu 

_oo 


g(u)  = f(t)e"1wtdt  . 

_oo 


This  is  accomplished  using  the  substitutions 


L-L’  +f 


F(L ' t ) = X ( L ' - J ) 
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(3c) 


G(2TrK/Nt0)e"11T^K"1^  = NtQC(K) 

and  by  changing  the  limits  on  K to  (1-  N/2,  N/2)  which  can  be  done 
because  C(K)  is  periodic,  C(K)  = C*(N-K+2).  These  substitutions  allow 
Equations  (1)  to  be  written  in  the  form 

| i jjr  (K-l)(L-l)t0 

F(uo> ' b * „ G(KZ  /Nto)e  0 Hr  (,a) 


f - (K-l)(L'-l)t0 

G(2tt K/Nt  ) = l F(L' t )e  0 t . (4b) 

0 , , , N 0 0 


The  substitutions  in  Equations  (3)  serve  two  purposes.  They  allow 
the  identfication  of  the  correspondence  between  the  discrete  and  con- 
tinuous variables,  as  in  the  following  expressions; 


t n,  Nt 

0 

(5a) 

dt  n-  t 

0 

(5b) 

w 'v  2ir«/Nt0 

(5c) 

doj  n.  2iT/Nt„ 
0 

(5d) 

f(t)  * F(L*t0)  = X(L * ) 

(5e) 

g(w)  ^ G(2ttK  /Nt0)e11T^K_1  ^ 5 Nt  C(K) . 

(5f) 

The  second  purpose  is  for  completeness.  It  is  to  shift  the  summation 
limits  to  make  them  symmetric.  This  is  done  in  the  frequency  domain 
by  changing  the  summation  limits  from  (1,N)  to  (l-N/2,  N/2)  (allowed 
since  the  representation  is  periodic)  and  by  translating  the  data 
string  to  make  it  symmetric. 


f 
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It  is  assumed  that  the  data  is  sampled  sufficiently  rapidly  so 
that  the  summations  accurately  represent  the  respective  integrations. 

I | 

The  autocorrelation  R(Lt0),  and  the  power  spectral  density, 

S(K2  /Nto)  are  defined  respectively  in  Equations  (6).  The  angular 
brackets  indicate  ensemble  averaging. 


R(Lt0)  = <F(It0)F*((I+L-l)t0)> 


1 

2* 


Y 

l-N/2 


i (K-1)(L-1) 

S(K2  /Nt0)e 


(6a) 


(6b) 


These  are  to  be  compared  with  the  corresponding  equations  for  the  con- 
tinuous case: 


R(trt2)  = <f(t1)f*(t2)> 


R(tl»t2)  = 2tt 


W-w)e 


iw(t-|-t2) 


We  also  have  for  the  continuous  case. 


(7a) 

(7b) 


<g(“1)g*(u2)>  = W(w-|)6(wi"u2^  • (7c) 

Clearly  the  spectral  density,  W(w)  corresponds  to  S(K2ir/Nt0)  in 
the  discrete  case. 

The  power  spectral  density,  S,  for  the  discrete  case  is  related 
to  the  transform  of  the  data,  C(K)  by  Equation  (8). 


-i*(K,-K?) 

Nt0<C(K12ir/Nt0)C*(K22Tr/Nt0)  >e 


S ( K2tt/N  tQ ) 


(8) 


Equation  (8)  is  the  discrete  counterpart  of  the  continuous  form  in 
Equation  (7c).  Thus  for  any  single  data  transform  we  take  for  com- 
putational purposes 


L 


6 


V 


2 


(9) 


S(K2u/Nt0)  = T|C(K2ir/Nt0) 


where 


T = Nt0  (10) 

We  note  that  there  are  only  N/2  unique  values  for  the  spectrum, 
S(k2  /Nt0) . This  follows  because  S(K2ir/Nt0)  is  symnetric  about  zero 
frequency.  Thus  the  frequency  dynamic  range  is  ir/t0,  covering  the 
range  (2ir/Nt0,ir/t0) . 

B.  Trend  Removal 

To  improve  our  estimates  of  the  very  low  frequency  components  of 
our  data,  it  is  beneficial  to  remove  the  DC  and  linear  trends  [4] 
before  calculating  the  power  spectrum. 

The  suposition  that  the  data  of  interest  are  modulated  by  a 
linear  trend  is  expressed  as; 


= Po  + 


Vi  + xi’ 


(lla) 


where  the  Xi  are  the  original  data,  B0  and  Bi  are  respectively  esti- 
mates of  the  DC  and  linear  trend,  and  the  x.j  are  the  data  of  interest. 
It  is  easily  shown  [1]  that  for  a minimum  mean  square  error  fit. 


P1  = 


EVxi 


Nit?  - (It.)2 


5o  ■ l <£* i - 51  h > 


(Hb) 


(11c) 


Without  loss  of  generality  we  can  assume  t.=(i-l).  In  this  case 
Equations  (11)  simplify  somewhat  to  yield. 
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8,  = [£(i-l)x.  - Zxi]  (12a) 

1 N(N2-1)  1 1 1 


N(N-1 ) 
2 


] • 


(12b) 


The  data  of  interest  are  then  given  by 


xi 


xi 


60  - Bl(i-1) 


(13) 


The  reason  for  subtracting  linear  trends  is  to  improve  low-frequency 
estimates  when  the  data  string  is  not  sufficiently  long  to  give  in- 
formation about  the  lowest  frequencies.  Thus  a linear  trend  could  be 
part  of  a sine  wave  with  period  longer  than  the  data  duration.  The 
effect  of  a trend  is  to  augment  the  values  of  the  lowest  few  harmonics 
by  adding  in  the  Fourier  series  spectrum  of  a linear  waveform  thus 
giving  incorrect  values.  The  existence  of  the  lower  frequencies  can 
be  indicated  by  the  shape  of  the  lower  frequency  portion  of  the  spec- 
trum. If  the  spectral  values  continue  to  rise  as  the  frequency  de- 
creases right  down  to  the  lowest  harmonic,  then  the  existence  of  even 
lower  frequencies  would  seem  quite  possible. 


C.  Hanning  and  Compensation  for  Non-Unity  Energy 


The  use  of  the  finite  discrete  Fourier  transform  pair  (Equations 
(3a, b))  implies  that  we  are  looking  at  a finite  section  of  an  infinite 
string  of  data  [7].  That  is,  we  are  viewing  our  data  through  a rec- 
tangular data  window  of  unity  height  and  length  T.  Multiplication  of 
our  data  by  this  window  results  in  the  true  spectrum  being  convolved 
with  a sin(x)/x  function*.  Since  this  function  has  significant 
sidelobes,  we  wish  to  somehow  alter  our  rectangular  window  so  that 
the  resulting  convolution  is  with  a function  with  significantly  lower 
side  lobes.  Such  a data  window  commonly  known  as  the  Hanning  window  [4] 
is  given  by 


cos 


(L-l-N/2)* 

N/2 


(14) 


*This  error  is  often  referred  to  as  leakage. 
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Equivalent  to  multiplying  our  data  by  this  window  and  transforming,  is 
simply  calculating  [2]  in  the  frequency  domain 


C(K)  H2  1 


1 

4 


- \ Re(C2) 

c + 1 c 

LK-1  2 LK 


K=1 


4 CK+1 ’ K_2,35‘ 


N 

2 


(15) 


Note  that  by  use  of  our  data  window  as  in  Equation  (14),  we  have  de- 
creased the  energy  in  our  original  data.  Therefore  to  obtain  an 
unbiased  estimate  of  the  true  power  spectrum  [3],  we  require  that  the 
Hanned  estimation  on  the  right  hand  side  of  Equation  (9)  be  multiplied 
by  8/3.  See  Appendix  A for  a demonstration  of  this  factor. 

D.  Pre-Whitening  and  Post-Darkening 

We  next  consider  the  problem  of  estimating  the  harmonics  of  a 
spectrum  which  has  a large  dynamic  range,  i.e.,  one  which  is  steeply 
sloping.  Recall  from  Section  C that  the  spectrum  we  calculate  is 
really  the  true  spectrum  convolved  with  the  transform  of  our  data 
window.  If  the  true  spectrum  is  steeply  sloping  even  though  the  side 
lobes  of  the  transform  of  the  data  window  are  small,  the  product  of 
the  harmonics  of  the  side  lobe  and  the  off-center  harmonic  (which  is 
much  greater  than  the  center  harmonic)  of  the  true  spectrum  may  be 
comparable  to  the  estimate  of  the  center  harmonic.  Therefore  to 
reduce  this  problem,  we  resort  to  the  technique  of  pre-whitening  [4] 
which  effectively  evens  out  the  true  spectrum.  Such  a pre-whitening 
scheme  is  implemented  by  calculating 


x 


i 


® Xj  i»  i— 2,3,***,N 
; i=i , 


(16) 


where  a is  the  pre-whitening  parameter.  Some  a priori  knowledge  of 
the  general  shape  and  dynamic  range  of  the  true  power  spectrum  is 
needed  to  choose  a value  for  a.  Use  of  expression  (16)  results  in 
the  true  power  spectrum  being  multiplied  [4]  by 


1-a  e 


i 2tt(L-1  ) 
N 


(17) 
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This  function  is  shown  plotted  in  Figure  3 for  several  representative 
values  of  a.  For  the  high  frequency  range,  a=  0.975  seemed  to  be  a 
good  choice  for  our  data. 

So  to  summarize,  our  procedure  for  power  spectrum  estimation 
amounts  to; 

1)  removing  DC  and  linear  trends  (Equations  (12)  and  (13)) 

2)  prewhitening  (expression  (16)) 

3)  calculating  the  voltage  spectrum  C(K)  (Equation  (3b)) 

4)  hanning  the  spectral  components  C(K)  (expression  (15)) 

5)  calculating  the  power  spectral  harmonics  S(K)  (Equation  (3a)) 

6)  compensating  for  non-unity  energy  in  the  data  window  by 
multiplying  the  S(K)  by  8/3 

7)  post-darkening  by  dividing  the  power  spectral  estimates 
S(K)  by  expression  (17). 

E.  Ensemble  Averaging,  Confidence  Intervals 

We  now  come  to  the  question  of  how  adequately  our  spectra  estimate 
the  true  power  spectrum.  If  we  assume  we  have  calculated  several 
spectra  (of  successive  data  strings  of  period  T)  then  assuming  that 
all  the  data  are  from  a random  stationary  process,  we  can  average  at 
each  spectral  harmonic  [5],  over  all  the  spectral  realizations.  That 
is  if  the  ith  harmonic  of  the  jth  spectral  realization  is  denoted  as 
S- .,  then  our  average  spectrum  will  be  given  by 

1 J 

■ r j,  sij  <>8> 

s)  ’ 

where  S.  denotes  the  estimate  of  the  ith  harmonic  of  the  true  spectrum 

C 


If  we  now  assume  that  the  original  data  were  samples  from  a 
gaussian  random  process  with  zero  mean  and  power  spectrum  Si  then  the 
statistic 


Y = 


(19) 


has  an  approximate  chi-square  distribution  [6]  with  2L  degrees  of 
freedom  (x^(2L)).  The  factor  of  2 arises  because  Si  is  the  modulus 
of  a complex  random  variable  (2  degrees  of  freedom). 
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HARMONIC  NO. 


Figure  3 


Pre-whitening  power  transfer  function  for  various 
values  of  the  pre-whitening  parameter  a. 


With  the  assumption  that  Y is  x^(2L)  we  can  calculate  the  proba- 
bility that  the  true  spectrum  lies  within  a certain  region  of  our 
estimate  Such  a probability  is  expressed  as 

P(a  < Y < b)  = 1 - Y fcOa! 


2L;  \ 


2L ; 1 - £ 


1 - Y 


(20b) 


where  x 


2L ; \ 


and  x 


2L;  1 - j 


are  respectively  the  y/2  and  1-y/2  per- 


centile points  of  a x^  distribution  with  2L  degrees  of  freedom.  With 
some  algebraic  manipulations,  Equation  (20b)  becomes 


± < ii<  1L 

2 - 2 
< Y S,  X Y 

2L;1-  J 2L;  J 


= 1-Y  • 


Equation  (21)  states  that  on  the  basis  of  the  available  data,  with 
probability  1-y,  the  true  spectral  value  Si  lies  within  the  limits 

% S.  and  \ S.  . (22) 

x Y x 

2L;1-  \ 2L;  \ 

Figure  4 is  a display  of  these  confidence  limits  plotted  vs  the  number 
of  spectra  in  the  average. 

This  concludes  our  discussion  of  several  classical  techniques 
employed  in  the  digital  estimation  of  power  spectra  for  the  high- 
frequency  range.  We  have  described  the  process  of  spectrum  calcu- 
lation via  the  FFT,  the  techniques  of  trend  removal,  banning  and 
compensation  for  non-unity  energy  in  the  hanning  window,  pre-whitening/ 
post-darkening,  and  a goodness  of  estimate  criterion  for  our  final 
power  spectrum. 

The  spectra  calculated  by  the  methods  outlined  in  this  section 
were  taken  to  be  bandlimited  to  w = ± 9.629  x 10^  rad/sec  and  have 
the  first  harmonic  at  Aw  = 47.017  rad/sec.  In  the  next  section  we 
will  give  a procedure  for  estimating  frequencies  lower  than  47  rad/sec. 
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EQUA- PROBABILITY  LIMITS 


Figure  4.  Normalized  confidence  limits  as  a function  of  the  number  of  ensemble  members. 


SECTION  IV 
LOW  FREQUENCY  RANGE 

We  now  present  a technique  for  estimating  frequency  components 
lower  than  the  first  harmonic  of  the  high  frequency  range. 

If  we  recall  Equation  (5d) 


(23) 


where  f = l/tQ 

we  see  that  there  are  two  ways  to  estimate  lower  frequencies:  (1)  use 

a larger  FFT,  i.e.,  make  N larger  or  (2)  decrease  f$.  Since  for  prac- 
tical purposes  we  are  limited  by  our  equipment  to  an  FFT  with  4096 
points,  we  must  decrease  the  sampling  frequency.  Merely  decreasing 
the  sampling  frequency,  however,  would  result  in  aliasing  of  the  es- 
timated spectrum.  Therefore  to  avoid  aliasing  we  must  first  apply  a 
low-pass  digital  filter  [4]  to  our  data  and  then  decimate  (take  every 
4th  or  every  10th  or  every  64th  etc.  point)  the  resulting  band-limited 
data. 


I 


Such  a smoothing  operation  (low-pass  filtering)  can  be  performed 
as 


Y 


n 


n+«.-l 

I 


J=n 


(24) 


where  * is  to  be  chosen  later.  It  is  easily  shown  that  application  of 
Equation  (24)  results  in  the  power  spectrum  being  multiplied  by 

'sin  MK-U/N  2 
sTn  ff(K-l )/N 

This  function  effectively  band-limits  the  spectrum  to  harmonics 
one  through  N/t  at  which  point  the  transfer  function  has  its  first  zero. 

We  may  also  reduce  the  side  lobes  of  this  transmission  function 
by  using  successive  smoothing  operations,  e.g.. 


where  p^t.  Obviously  these  two  successive  smoothing  operations  multiply 
the  true  spectrum  by 
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J 


2 


(25) 


sin  iu(K-l  )/N  1^| 

sin  prr ( K-l  )/N 

sin  n ( K- 1 )/N 

sin  tt ( K-l  )/N 

With  the  smoothing  operation  thus  accomplished  we  can  then  decimate 
our  data  without  folding  harmonics  N/t  through  N/2  - 1 onto  harmonics 
zero  through  N/£  - 1. 


Then  one  can  go  on  with  the  general  procedure  outlined  previously. 
After  performing  these  calculations,  the  resulting  spectra  are  then 
simply  divided  by  expression  (25). 


sin  £tt(k-1  )/(nN) 
sin  tt ( K- 1 )/(nN) 


sin  ptt ( 

;k-ii 

i/i 

[nNl 

1 

sin  tt | 

!k-i] 

i/i 

1 

(26) 


where  n is  the  decimation  factor.  This  division  operation  deserves 
one  word  of  caution:  Division  by  the  small  magnitude  of  the  filtering 

transfer  function  near  the  high  frequency  end  of  the  pass-band  em- 
phasizes the  effects  of  noise.  Therefore  these  last  few  high  frequency 
estimates  should  be  discounted. 

For  the  spectra  calculated  herein  we  chose  to  smooth  successively 
by  groups  of  7 and  8.  The  power  spectrum  of  these  smoothing  operators 
is  shown  in  Figure  5.  With  the  smoothing  thus  accomplished,  we  then 
decimated  by  4,  i.e.,  we  kept  only  1/4  the  number  of  points.  Figure 
5 also  shows  the  resulting  aliasing.  Going  back  to  Equation  (4)  of 
the  previous  section  we  see  that  we  have  decreased  the  harmonic  spacing 
(and  the  bandwidth)  by  a factor  of  four. 

The  advantage  of  this  smoothing  and  decimation  technique  is  that 
one  can  repeat  the  operation  many  times.  Thus  we  can  estimate  ar- 
bitrarily low  (limited  only  by  the  total  amount  of  data)  frequencies. 

Smoothing  techniques  as  evidenced  by  Equation  (24)  are  of  a general 
class  known  as  nonrecursive  filtering.  The  more  elegant  recursive  fil- 
tering employs  a linear  combination  of  the  form 


Y 

n 


1=0  a,V'  * 1=1  b,V"-1 


With  the  proper  choice  of  weights  a-j  and  b-j  this  summation  yields  the 
digital  analog  of  more  sophisticated  filters,  e.g.,  Butterworth, 
Chebyshev.  Reference  10  contains  a comprehensive  discussion  of  these 
techniques. 


15 


This  ends  the  theoretical  background  for  the  estimation  of  micro- 
thermal  fluctuation  power  spectra.  In  the  next  section  we  consider  the 
selection  of  spectral  ranges,  i.e.,  data  record  lengths  and  decimation 
factors  for  the  desired  confidence  intervals. 


O 

O 


Figure  5.  Power  transfer  function  for  smoothing 
by  groups  of  7 and  8. 
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SECTION  V 


SELECTION  OF  SPECTRAL  RANGES 

We  now  consider  a question  not  mentioned  previously.  That  is, 
given  a very  large  amount  of  data  how  should  it  be  divided  in  order 
to  obtain  the  best  spectra  the  most  efficiently?  This  question  is 
faced  every  time  one  analyzes  a data  sequence  of  more  points  than  the 
largest  available  FFT.  The  answer  determines  the  precision  of  the 
spectral  estimates  and  the  required  computer  time  and  must  certainly 
be  determined  before  starting  the  procedures  of  Sections  III  and  IV. 
Here  we  present  recommendations  to  assist  in  the  solution  of  this 
problem  for  any  particular  case. 

The  result  of  the  procedures  described  in  this  chapter  is  a 
composite  spectrum  covering  several  different  ranges.  These  will  be 
obtained  from  data  sequences  of  different  time  durations.  The  higher 
the  frequency,  the  shorter  the  time  duration  of  data  which  is  required. 
The  data  sequence  for  each  spectral  range  is  divided  into  sections, 
the  spectrum  is  obtained  for  each,  and  all  the  spectra  are  averaged 
point  by  point  to  provide  a statistically  reliable  estimate.  For 
all  but  the  highest  frequency  range  the  data  sequence  has  been  formed 
from  the  filtering  and  decimation  procedures.  The  lowest  frequency 
ranges  will  have  no  statistical  averaging  since  the  data  have  been 
decimated  to  the  point  where  there  are  sufficient  points  for  only  a 
single  FFT.  As  a rule  there  is  reduced  statistical  reliability  at 
the  low  frequency  end  of  the  spectrum. 

We  assume  as  starting  information  that  the  data  sequence  has  a 
total  time  duration  T],  is  bandlimited  and  has  been  digitized  at  a 
rate  f0,  which  is  greater  than  or  equal  to  the  Nyquist  rate.  We 
expect  as  a final  result  a composite  power  spectrum  with  a lowest 
frequency  around  1/Ti  and  a highest  frequency  of  f0/2.  Four  points 
should  be  considered  in  the  section  of  spectral  divisions: 

1.  the  duration  of  the  data  sections  for  each  spectral  range 

2.  the  number  of  spectra  in  the  ensemble  average  for  each 
spectral  range 

3.  the  factor  by  which  to  decimate 

4.  the  treatment  of  the  lowest  frequency  range. 

The  question  of  the  duration  required  for  the  data  sections  will 
be  considered  first.  The  discussion  is  aimed  primarily  at  the  high 
frequency  range,  however  it  is  applicable  with  the  proper  frequency 
scaling  to  all  frequency  ranges.  The  basic  issue  is  whether  to  employ 
one  large  or  several  small  FFT's.  In  the  former  case  there  would  be 
fewer  spectral  ranges  covering  the  same  span  of  frequencies.  The 
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answer  depends  upon  the  computer  time  required  to  perform  large  FFT's 
and  a few  filtering  and  decimation  procedures  as  opposed  to  small  FFT’s 
and  many  filtering  and  decimation  procedures.  Since  for  many  computers 
(as  was  the  case  here)  the  filtering  and  decimation  procedure  requires 
significant  amounts  of  time,  choosing  the  largest  possible  FFT  gives 
the  more  efficient  procedure.  If  the  largest  FFT  required  N points 
then  the  highest  frequency  range  will  span  the  frequencies  f0/N  to 
f0/2  as  shown  in  Figure  6. 

Now  consider  the  question  of  statistical  reliability.  The  desired 
precision  is  used  in  conjunction  with  Figure  4 which  gives  the  ratio 
of  confidence  interval  (normalized  to  the  mean  value)  as  a function 
of  M,  the  number  of  spectra  averaged.  In  addition  it  simplifies  the 
decimation  operation  if  M is  chosen  to  be  a power  of  2.  For  our  case 
M=32  seemed  a good  choice. 

The  next  issue  is  the  decimation  rate.  The  guideline  here  is  to 
decimate  by  a sufficient  factor  so  that  there  is  only  a small  overlap 
between  the  adjacent  spectral  ranges.  An  overlap  of  a factor  of  8 
or  16  seems  to  be  reasonable.  The  factor  of  8 gives  at  most  eight 
points  for  comparison. 

For  decimation  by  a factor  of  2n  the  next  to  highest  frequency 
range  will  run  between  f0/N2n  and  fo/2(2)n,  and  after  decimating  m 
times  the  resulting  frequency  range  will  be  between  f0/N2nm  and 
f0/2(2)nm.  The  next  to  highest  frequency  range  will  require  a data 
duration  of  MN(2n)/f0  seconds  and  after  decimating  m times  the  re- 
quired data  duration  will  be  MN(2nm)/f0  seconds.  Finally  we  require 
that  after  m decimations  each  by  a factor  of  2n  we  are  left  with 
exactly  one  record  of  N points. 

The  last  point  is  the  treatment  of  the  lowest  frequency  range. 

This  is  important  because  in  the  low  frequency  ranges  the  duration  of 
the  individual  sections  is  sufficiently  long  to  reduce  the  number  of 
sections,  thus  giving  reduced  statistical  reliability.  Indeed  in 
the  last  range  there  will  be  only  one  section  of  duration  T] . The 
single  spectrum  yields  minimum  statistical  confidence.  For  this  fre- 
quency range  confidence  intervals  may  be  improved  by  dividing  the  N 
remaining  data  points  into  several  smaller  sub-sections  and  performing 
a smaller  FFT  on  each.  This  increases  the  lowest  estimable  fre- 
quency by  a factor  equal  to  the  number  of  subsections.  Therefore 
we  see  that  as  one  attempts  to  estimate  lower  and  lower  frequencies 
(assuming  the  data  is  of  limited  extent)  the  confidence  intervals  must 
necessarily  increase. 

The  low  frequency  problem  can  be  minimized  in  many  cases  by 
choosing  the  digitizing  frequency  f0  so  that  there  are  exactly  N data 
points  after  the  final  filtering/decimation  operation  and  exactly  M 
sections  of  N points  each  before  the  final  smoothing  and  decimation 
operation.  This  requires  that  M,  the  number  of  sections  averaged  be 
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Figure  6.  Frequency  ranges  for  various  decimation  factors. 


equal  to  2n,  the  decimation  factor,  and  that  the  digitization  frequency 
be  f0=NK/T-|,  where  k=2nm  is  the  original  number  of  sequences  of  N points 
each  in  the  unfiltered  data. 

22 

As  ai  example  suppose  we  have  NK=2  data  points  and  the  largest 
available  FFT  has  N=2 ' 2 points.  Then  we  have  K=2^ 0=1024  sequences  of 
4096  points  each.  The  desired  confidence  intervals  require  averaging 
over  M=2^=32  spectra.  Consequently  we  must  have  K=2nm  0r  n=Log2M=5 
so  that  m=2.  That  is  we  need  m=2  decimations  eacn  by  a factor  of 
2n=32.  For  the  highest  frequency  range  we  need  only  calculate  32  of 
the  possible  1024  spectra  to  attain  the  desired  precision.  After  the 
first  decimation  by  2$  we  are  left  with  32  records  of  N points  each, 
and  after  the  second  decimation  with  exactly  one  record  of  N points. 

This  example  gives  an  overlap  between  successive  spectral  ranges  of 
2&=64  points. 

This  completes  the  discussion  of  the  selection  of  spectral  ranges. 
The  procedure  outlined  herein  offers  a straightforward  approach  aimed 
at  providing  the  desired  precision  of  the  spectral  estimates  con- 
sistent with  minimum  computer  time. 


SECTION  VI 
RESULTS 


We  now  present  the  results  of  the  procedures  which  were  outlined 
in  Sections  III,  IV,  and  V. 

Our  results  consist  of  several  plots,  each  marked  with  a group 
number,  a range  number  and  an  average  over  K files.  The  group  number 
refers  to  a particular  8.36  second  segment  of  data  (of  which  there  are 
15).  The  range  number  n refers  to  the  spectra  calculated  from  data 
which  have  been  decimated  by  a factor  of  4(n'N.  Each  file  refers  to 
a single  4096  point  FFT,  and  the  number  K denotes  the  number  of  spectra 
averaged  to  produce  the  plot  (K=4(n"')). 

Plots  7 and  8 show  the  spectra  calculated  for  ranges  1 and  2 
respectively.  We  note  the  presence  of  the  noise  spikes.  The  first 
spike  is  at  approximately  400  Hz,  the  second  and  third  spikes  are 
respectively  3rd  and  5th  harmonics  of  400  Hz.  The  remaining  spikes 
are  all  approximately  harmonically  related  and  as  such  can  be  dis- 
counted as  legitimate  data. 

Plots  9 through  12  show  the  complete  results  for  group  1.  The  noise 
spikes  evident  in  ranges  1 and  2 have  been  removed.  Note  that  in  the 
regions  of  overlap  between  successive  ranges  there  is  very  good  agree- 
ment between  spectral  estimates. 

Next,  plots  13  through  15  are  the  same  as  plots  three  through  five 
with  the  exception  that  we  have  also  plotted  the  90%  confidence  in- 
tervals. 

Finally,  Figures  16  are  composites  of  the  spectra  from  all  fre- 
quency ranges  averaged  over  the  full  2.14  minutes.  As  a check  on  the 
spectrum  calculated  herein,  the  signal  variance  over  the  entire  2.14 
minute  run  was  estimated  graphically  from  the  spectrum  of  Figures  16-a 
and  compared  with  the  variance  calculated  from  the  raw  data.  The  two 
estimates  were  found  to  differ  by  less  than  2%. 

By  averaging  spectra  over  the  entire  stretch  of  data  we  have 
obviously  made  the  assumption  that  the  process  giving  rise  to  the 
temperature  fluctuations  is  wide-sense  stationary.  That  this  might 
be  the  case  is  indicated  by  the  map  shown  in  Figure  17. 

The  position  and  flight  length  during  which  the  data  were  taken 
is  indicated  by  the  large  black  arrow.  As  can  be  seen,  there  are  no 
large-scale  geographical  features  in  the  area  of  data  acquisition 
which  would  destroy  the  stationarity  of  the  data.  The  map  also  shows 
surface  features  which  could  have  a profound  effect.  Therefore  this 
point  should  be  checked  on  each  data  run  until  more  is  known  about  the 
effects  of  terrain  features. 
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Figure  18  shows  the  spectrum  obtained  with  the  standard  procedure 
used  at  Kirtland  Air  Force  Base.  The  curve  has  been  normalized  so 
that  the  area  under  the  curve  is  equal  to  2ira2  where  a?  is  the  variance 
estimate  obtained  in  this  report.  This  normalization  corrects  for  a 
discrepancy  in  spectral  magnitudes.  Comparison  of  Figures  16-a  and  18 
shows  that  the  spectra  agree  quite  closely. 
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Figure  9.  Same  plot  as  in  Figure  7 with  frequency 
spikes  removed. 


SPECTRUM  (DE6*k2  SEC) 


RANGE  N0.  2 


AVERAGE  OVER  16  FILES 


GROUP  NO.  1 


-3  ' -2  1 -I 

10  10  10 


Figure  10.  Sa 
sp 


Figure  11.  Third  frequency  range  spectrum. 
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Figure  13.  Spectrum  with  90%  confidence  interval. 
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Figure  14.  Spectrum  with  90%  confidence  interval. 
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Figure  16-b.  Spatial  spectrum  corresponding  to  Figure  15-a 
(velocity  = 155  m/sec.). 
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Figure  17.  The  arrow  at  the  middle  right  of  the  map  indicates  the  direction 
and  distance  flown  during  data  acquisition. 


SPECTRUM 


Figure  18.  Comparison  of  spectrum  calculated  at  Kirtland  Air  Force 
Base  (solid  lines)  and  the  spectrum  calculated  herein 
(points) . 
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SECTION  VII 

SUMMARY  AND  CONCLUSIONS 


This  report  has  considered  the  calculation  of  temporal  spectra 
of  microthermal  fluctuations  measured  aboard  a jet  aircraft.  The 
calculations  were  performed  as  a check  on  spectra  obtained  at  Kirtland 
Air  Force  Base. 

Section  II  discussed  the  digitization  of  the  analog  data  plus 
some  preliminary  processing  such  as  removal  of  obviously  erroneous 
data  spikes. 

Section  III  covered  the  various  techniques  employed  in  calculating 
the  power  spectrum  by  use  of  the  F.F.T.  algorithm.  Such  techniques 
include  trend  removal,  Hanning  and  compensation  for  non-unity  energy, 
prewhitening  and  post  darkening  and  calculation  of  confidence  in- 
tervals. Also  contained  in  this  section  is  a discussion  of  the  analogy 
between  the  equations  implemented  by  the  F.F.T.  and  the  integral 
Fourier  transform  pair. 

The  calculation  of  spectra  for  arbitrarily  low  frequencies  is 
demonstrated  in  Section  IV.  In  addition  to  the  technique  outlined 
above,  use  is  made  herein  of  digital  filtering. 

Selection  of  sampling  rates,  and  spectral  ranges  which  will  yield 
the  desired  confidence  levels  are  discussed  in  Section  V. 

The  results  of  the  calculations  are  contained  in  Section  VI  in 
the  form  of  several  graphs. 

As  a result  of  the  present  analysis  several  heretofore  unknown 
features  of  the  data  were  uncovered  and  some  suggestions  for  future 
data  analysis  were  made. 

It  can  be  concluded  that  there  was  good  agreement  in  shape  and 
close  but  not  exact  agreement  in  magnitude  between  the  spectra  calcu- 
lated herein  and  those  calculated  at  Kirtland  Air  Force  Base. 
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APPENDIX 


In  this  appendix  we  demonstrate  the  decrease  in  relative  energy 
experienced  when  a non-rectangular  data  window  is  employed.  For  the 
Hanning  window  the  decrease  is  shown  to  be  by  a factor  of  3/8. 

The  spectrum  as  expressed  by  Equation  (lb)  in  the  text  is 

, N -i  (L-1)(K-1) 

C(K)  = l l X(L)e 
N L=1 

The  total  energy  is  thus  given  by 


K=1 


I C ( K)  | 2 = ^ l l l X(L)e 
K=1  L=1  M=1 


-i  (L-1)(K-1)  1 (L-1)(M-1) 

N X*(M)e  N 


! N N N -1  (L-1)(K-M) 

K l l X(L)X*(M)  l e N 


N L=1  M=1 


K=1 


N N 


= K l l X(L)X*(M)[N6  M] 
L=1  M=1 


Then 


N o i N , 

E = l |C(K)|2  = Jr  l |X(L)|2. 
K=1  n L=1 


For  the  rectangular  data  window  the  relative  energy  is  (X(L)=1); 
L=1 ,2, • • • ,N) 


is 


N 


er=H  i'-l. 


L=1 


Now  the  Hanning  window  as  expressed  by  Equation  (14)  in  the  text 


X ( L ) 


*[ 


! + cos 

i + cos  N/2 
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And  the  relative  energy  is 


E = — 
lH  4N 


J [,  + cos  (H-H/j)' 

L i * cos  N/2 


L=1  L 


J_  ? f,  + 2 cos  (l-l-i./?)"  * cos2  (L-T-N/2)^ 
4N  , , L!  c C0S  N/2  C0S  N/2  J 


, ' , (f,'f  U-l-K/2)^-,  p.  (L-l-N/2) 

4N  L-l  "»  lU  2 


J_  y e 


i ^ (L-l-N/2)  -i  ^ (L-l-N/2) 
n +2  + e 


+ l ]_  ? Jr  <L-" 

4 ' 4N  J,  e 


1 y -'r(L-i) 

® j, e 


i_  y j r ,L-,)  i j_  y 

™ i e 8 16N  e 


+ tcm  1 e 
16N  L=1 


E .1+1-1 

lH  4 8 8 


Q.E.D. 
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METRIC  SYSTEM 


BASE  UNITS: 

Quantity 


Unit 


length 

mass 

time 

electric  current 
thermodynamic  temperature 
amount  of  subatance 
luminous  intensity 

SUPPLEMENTARY  UNITS: 

plane  angle 
solid  angle 

DERIVED  UNITS: 
Acceleration 

activity  (of  a radioactive  source) 

angular  acceleration 

angular  velocity 

area 

density 

electric  capacitance 

electrical  conductance 

electric  field  strength 

electric  inductance 

electric  potential  difference 

electric  resistance 

electromotive  force 

energy 

entropy 

force 

frequency 

illuminance 

luminance 

luminous  flux 

magnetic  field  strength 

magnetic  flux 

magnetic  flux  density 

magnetomotive  force 

power 

pressure 

quantity  of  electricity 
quantity  of  heat 
radiant  intensity 
specific  heat 
stress 

thermal  conductivity 
velocity 

viscosity,  dynamic 

viscosity,  kinematic 

voltage 

volume 

wavenumber 

work 


SI  PREFIXES: 


metre 

kilogram 

second 

ampere 

kelvin 

mole 

candela 


radian 

steradian 


metre  per  second  squared 

disintegration  per  second 

radian  per  second  squared 

radian  per  second 

square  metre 

kilogram  per  cubic  metre 

farad 

siemens 

volt  per  metre 

henry 

volt 

ohm 

volt 

joule 

loule  per  kelvin 

newton 

hertz 

lux 

candela  per  square  metre 
lumen 

ampere  per  metre 

weber 

tesla 

ampere 

watt 

past  al 

coulomb 

joule 

watt  per  steradian 

loule  per  kilogram-kelvin 

pascal 

watt  per  metre-keivin 
metre  per  second 
pascal-second 
square  metre  per  second 
volt 

cubic  metre 
reciprocal  metre 
joule 


Multiplication  Factors 


1 000  000  000  000 

* 

10” 

1 000  000  000 

■ 

10’ 

1 000  000 

C 

10* 

1 000 

e 

10’ 

100 

* 

10> 

to 

* 

10' 

0 1 

= 

10-' 

0 01 

a 

irr' 

0 001 

* 

10*> 

0 000  001 

s r 

10  * 

0 (Kill  000  (XU 

a 

10- * 

o (too  ooo  (xxi  ooi 

* 

10-  ” 

o ooo  ooo  ooo  txx)  ooi 

in-” 

0 000  000  000  000  (XX)  001 

* 

in'* 

* To  be  avoided  where  poaaible 


SI  Symbol  _ Formula 

m 

kg 

s 

A 


mol 

cd 

rad 

sr 

m/s 

(disintegration  )/s 

radis 

rad/s 

m 

kg/m 

F 

As/V 

S 

A/V 

V'm 

H 

Vs/A 

V 

W/A 
V/ A 

V 

W/A 

1 

N-m 

)/K 

N 

kg-m's 

Hz 

(cyclej/s 

lx 

Im/m 

cd/m 

Im 

cd-sr 

A/m 

Wb 

V-s 

T 

Wb/m 

A 

W 

l'» 

Pa 

N/m 

C 

A-s 

1 

N-m 

W'tr 

I'kg-K 

Pa 

N/m 

W/m-K 

m/s 

Pa-s 

m/s 

V 

W/A 

m 

|wsve)/m 

1 

N-m 

Prefix 

SI  Syn 

lera 

T 

R'X* 

C 

mega 

M 

kilo 

k 

hecto’ 

h 

deka* 

da 

decl* 

d 

i entl* 

c 

milli 

m 

micro 

M 

nano 

n 

pico 

(emln 

alto 

r 

■ 
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